# Alexander F. Gazmararian
# afg2@princeton.edu

# Load packages
library(tidyverse)
library(margins)
library(modelsummary)
library(emmeans)
library(sandwich)
library(estimatr)

# Packages for structural topic model
library(stm)
library(geometry)
library(Rtsne)
library(rsvd)

# Load functions
source("code/fun/book_theme.r")
source("code/fun/savefig.r")
source("code/fun/fix_txt.r")

# Set up coefficient names
source("code/fun/coefnames4tables.r")
coefnames <-
  c(
    coefnames,
    "lock_second" = "Lock-in order",
    "delegate_lead_treat"="Community delegation",
    "delegate_voice_treat"="Public input",
    "delegate_fund_treat"="Long-term funding",
    "delegate_lead_treat:delegate_voice_treat" = "Community delegation x Public input",
    "delegate_lead_treat:delegate_fund_treat" = "Community delegation x Long-term funding",
    "delegate_voice_treat:delegate_fund_treat" = "Public input x Long-term funding",
    "delegate_lead_treat:delegate_voice_treat:delegate_fund_treat" = "Community delegation x Public input x Long-term funding",
    "lockin_treat_poolGeneral Benefits" = "General benefits treatment",
    "lockin_treat_poolSpecific Benefits" = "Specific benefits treatment",
    "consensus_treat" = "Consensus treatment",
    "consensus_treat:PartySummaryNeither" = "Consensus treatment x Independent",
    "consensus_treat:PartySummaryRepublican" = "Consensus treatment x Republican"
  )

# Load data
g <- readRDS("data/NatQual_Winter22.rds")

# Specify baseline model
f.base <- y ~ age + Female + Black + Hispanic +  PartySummary + CollegeDegree + employfull + Rural_bin + Investment + ClimateImpact_tri + income5 + trust_index

# Estimate model
coalition <- list()
coalition[["Binary"]] <- lm(update(f.base, SecondReverse_bin ~ consensus_treat + .), g)
coalition[["Scale"]] <- lm(update(f.base, SecondReverse_scale ~ consensus_treat + .), g)
coalition[["Party"]] <- lm(update(f.base, SecondReverse_scale ~ consensus_treat * PartySummary + .), g)

## Create Figure 5.7 ----
m1 <- lm(reverse_likely ~ consensus_treat * PartySummary, g)
m1.mean <-emmeans(m1, ~ consensus_treat * PartySummary)
m1.mean <- data.frame(m1.mean)
m1.mean$consensus_treat <- ifelse(m1.mean$consensus_treat == 0, "Hidden", "Revealed")
m1.mean$lb90 <- with(m1.mean, SE * qt(.05, df) + emmean)
m1.mean$ub90 <- with(m1.mean, SE * qt(.95, df) + emmean)
m1.mean$PartySummary <- with(m1.mean, ifelse(PartySummary == "Neither", "Independent", as.character(PartySummary)))
m1.mean <- subset(m1.mean, PartySummary!="Independent")
m1.mean
plot.consensus <-
  m1.mean %>%
  ggplot(aes(x=consensus_treat,y=emmean,color=consensus_treat,label=scales::percent(emmean, accuracy=1))) +
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width = 0) +
  geom_errorbar(aes(ymin=lb90, ymax=ub90), width = 0, size = 1.75) +
  geom_point(size = 4) +
  facet_wrap(~PartySummary, scales = "free_x") +
  geom_text(hjust=-.75, size = 3, color = "black") +
  scale_color_grey() +
  book_theme +
  scale_y_continuous(label = scales::percent, limits = c(.49, .9), expand=c(0,0)) +
  labs(x="National Consensus for Compensation and Investments", y = "Likely reversal") +
  theme(legend.position="")
plot.consensus
savefig(plot.consensus, "5.7_figure_consensus", height = 2, filepath = "figures/")

## Create online appendix table ----
file <- "tables/ch5/ols_consensus.txt"
modelsummary(
  coalition,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefnames,
  gof_map = c("nobs", "adj.r.squared"),
  escape=FALSE,
  output="latex"
) %>%
  cat(., file = file)
fix_txt(file)

